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Abstract — This article concerns over-resolution methods for image reconstruction from data provided by a 
family of scanning instruments like the Herschel observatory. The work is centered on building a model of the 
instrument that faithfully reflects the physical reality, accurately taking the acquisition process into account so 
as to explain the data in a reliable manner. The inversion, i.e. the image reconstruction process, is based on a 
linear approach resulting from a quadratic regularized criterion and numerical optimization tools. The inversion 
process produces maps that are over-resolved with respect to the nominal resolution of the instrument. The 
application concerns the reconstruction of maps for the SPIRE instrument of the Herschel observatory. The 
numerical evaluation uses simulated and real data to compare the standard tool (coaddition) and the proposed 
method. We demonstrate a capacity to restore spatial frequencies over a bandwidth four times that possible with 
coaddition and thus to correctly show details invisible on standard maps. The proposed approach is also appUed 
to real data with significant improvement in spatial resolution. 

Key words. Techniques: image processing, data acquisition modelling, inverse problem, deconvolution, over- 
resolution, regularization, image processing. Methods: statistical, numerical. Astronomical instrumentation, 
methods and techniques 

1. Introduction 

Map making is a critical step in the processing of astronomical data from various imaging instruments (inter- 
ferometers, telescopes, spectro-imager. . .). A very large number of papers have been devoted to this question in 
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various literature and two recent special issues have been published |Leshem et al. 2010 2008| in the signal-image 
processing literature. Since the observed sky may contain structures of various scales, from extended emission to 
point sources, the challenge is to design reconstruction method delivering maps that are photometrically valid for 
the broadest range of spatial frequencies. 

For long-wavelength instruments, be they ground based (SCUBA/JCMT, LABOCA/APEX,. . . ), on-board bal- 
loons (Archeops, BLAST,. . . ) or space borne (IRAS, ISO, Spitzer, WMAP, Planck, Herschel,. . . ), the task is espe- 
cially challenging for two reasons. First, the resolution is poor at such wavelength. Second, the number of detectors 
in the focal plane of such instruments is limited, so, it is generally not fully sampled. Therefore, specific scanning 
strategies have to be defined, depending on the detector positions, and closely combined with a well designed 
image reconstruction method. 

The Herschel Space Observatory [Pilbratt et al. 2010| was launched in May 2009 together with the Planck 



satellite. It continuously covers the 55-672 iim spectral range with its very high resolution spectrometer HIFI | de 



|Graauw et al.||2010) and its two photometers / medium resolution spectrometers PACS | |Poglitsch et al.|[2"010) and 
SPIRE I Griffin et al. 2010| . With a 3.5 m primary mirror, Herschel is the largest space telescope launched to date. 
In order to take full advantage of the size of the telescope, the accurate representation and processing of the highest 
spatial frequencies presents a particular challenge. To this end, two step-by-step photometer pipelines have been 
developed by the instrument consortia by [Griffin et"aL]p008) for SPIRE and by |Wieprecht et al.| | |2009| for PACS: 
they produce flux density timelines corrected for various effects, calibrated and associated with sky coordinates 
(Level-1 products), then produce maps (Level-2 products). An important step is the correction of the 1// noise 
components, which can be correlated or uncorrected between bolometers. For SPIRE, a significant fraction of 
the correlated component is processed using the signals delivered by blind bolometers. For PACS, it is currently 
processed using different kinds of filtering. The glitches due to the deposit of thermal energy by ionizing cosmic 
radiation are flagged or corrected. Finally, the output can be simply coadded to produce "naive maps". Maximum 
likehhood approaches, namely MADmap jCantalupo et al.| [20T0) and SANEPIC JPatanchon et al.|[2008) , have 
also been developed to compute maps, using the spatial redundancy to correct for the 1 // noise. 

There are several drawbacks to these pipelines. First, as they work on a step-by-step basis, the global process 
suffers from the limitation of the worst step. Second, the information available in one step is not fully exploited 
by the others since only the final result of one step is handed over to the following. Third, the instruments and the 
telescope properties (mainly the diffraction) are not taken into account, so, the maps are unavoidably smoothed by 
the PSF, whereas the scanning strategy allows higher spatial frequencies to be indirectly observed. 



In order to overcome these limitations, we resort to an inverse problem approach | Idier 2008) that is founded 
on an instrument model and an inverson method. 



- It requires an instrument model that faithfully reflects the physical reality so as to distinguish, in the observa- 
tions, between what is caused by the instrument and what is due to the actual sky. To this end, an important 
contribution of our paper is an analytical model based on a physical description of the phenomena, as functions 
of continuous variables. Moreover, it includes scanning strategy, mirror, wavelength filter, feedhorns, bolome- 
ters and read-out electronics. As far as the resolution is concerned, the point is the following. On the one hand, 
the field of view is covered by hexagonally packed feedhom-coupled bolometers, and the distance between 
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adjacent bolometers is twice the PSF width, resulting in spectral aliasing. On the other hand, the scanning 
strategy introduces spatial redundancies contributing to a higher equivalent sampling frequency. So, it is cru- 
cial to properly take into account the scanning strategy and the whole instrument including aliasing in order to 
obtain over-resolution and reverse aliasing (see also the analysis of | Orieux et al. 2009| ). To the best of our 
knowledge, such a model has never been used in a map making method. 
- The inversion step constitutes an ill-posed problem |Idier[ |2008| due to the deficit of available information, 
and the ill-posedness becomes all the more marked as resolution requirement increases. The inversion meth- 
ods must therefore include other information in order to compensate for the deficits in the observations. Each 
reconstruction method is thus specialized for a certain class of maps (point sources, diffuse emission, superpo- 
sition of the two. . . ) according to the information taken into account. From this standpoint, the present paper is 
essentially devoted to extended emission. 

From the methodological point of view, the used inversion comes within the framework linear regularized 
method JTikhonov & Arsenin[ [ig??! [Andrews & Huntl|1977[ . It relies on a quadratic criterion involving an 
adequation measure (observed data vs model output) and a spatial smoothness measure. From a practical stand- 
point, we resort to a gradient based optimisation algorithm | Nocedal & Wright 2000| to compute the map. 
Moreover, in as much as it relies on two sources of information, the method is based on a trade-off tuned 
by means of an hyperparameter It is empirically set in the present paper and work in progress, inspired 
from I Orieux et al. 2010| , is devoted to the question of the hyperparameter and instrument parameters auto- 
calibration (myopic and unsupervised inversion). 



One of the most striking results of the present paper is the correct restoration of small-scale structures (wide 
band) whereas they are not detectable on naive maps. Such a result is reached thanks to the developped instrument 
model together with the used inversion: they jointly enable the proposed method to reduce instrument effects, over- 
take instrument limitations and restore high spatial frequencies. The latter point is closely related to the capability 
to reverse aliasing in an inversion process. 

In the image processing community, such capabilities are referred to as super-resolution [Park et aL] |2003[ . 
A lot of work has been devoted to this topic and the proposed paper is partly inspired by recent developments in 
this field. They are usually based on various (scene or camera) motion or scanning strategy. Some of them account 
for possible rotation plad & Feuer||1999) and/or a magnifying factor | |Rochefort et al.|[2006) . Other approaches 
introduce edge-preserving prior | Nguyen et al. 200 1[ Woods et al. 2006). Theses works rely on the description of 



the unknown object as a function of continuous variables, that is decomposed on pixel indicator basis | Hardie et al. 



1997 Patti et al. 1997| , on a truncated discrete Fourier basis | Vandewalle et al. 2007 1 , on a family of regularly 



shifted Gaussian functions |Rodet et al.||2008") , or spline family | |Rochefort et al.[ [20061 . Other approaches have 
been proposed, based on shift-and-add step | Farsiu et ar]|2004| followed by a deconvolution step | MoUna & Ripley 



1989). Finally, several contributions are devoted to the performance of super-resolution approaches |Champagnat 



[etaL] [20091 [Orieux et al.j [20091 • 

This paper is organized as follows. The instrument model describing the relationship between the measured 
data and the unknown sky is presented in Section[2] Section[3]details the method which is proposed to inverse the 
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data and compute high resolution maps. Finally, Section]?] presents experimental results, first on simulated data 
(Section |4~T] l then on real data (Section |4.2[ ). 



2. Instrument model 

The prime objective of the instrument model is to reproduce the data for a fixed sky by describing the physics of 
the acquisition. The sky, X{a, [3, A), is characterized by two spatial dimensions (a, f3) and one spectral dimension 
A. The closer the description of the acquisition is to reality, the better the inversion method is able to distinguish 
between information coming from the object and changes introduced by the instrument. 

The principle of faithfully representing the real world is difficult to conciliate with the practical constraints of 
numerical computing. In general, the reconstruction algorithms call on the instrument model many times in order 



to reproduce a data set resulting from a current sky. This is one of the differences with a simulator [Sibthorpe 



et al. 2009 Sibthorpe & Griffin 2006| which is designed to be run once per data set and is not included in a 



reconstruction algorithm. It is thus necessary to make choices or approximations to reduce computation time. 

2.1. Physical models 

2.1 .1 . Mode of observation 

To model telescope translations, we use a frame of reference defined by the instrument. The map present at the 
input is time dependent and can be written 

X{a, (3, \t)=X{a- p^{t), /? - pp(t), A) , (1) 

where a and (5 define the central angular position of the observation and {pa{t),Pi3{t)) the translations in the two 
directions as a function of time t. 

Here, we present only the Large map protocol for wide field of view. Data are acquired over a complete 
observation sequence composed of two practically perpendicular directions and several scans in one sense and the 
other for each of the two directions. The pointing acceleration and deceleration phases are outside the zone of 
interest and there is no rotation during the observation sequence. The pointing functions are thus written 

Pa{t)^Vat + Ca and pp{t) = Vpt + Cfi (2) 

for scanning at a constant velocity (wq, f^g). The pointing accuracy is of the order of a few seconds of arc. This 
protocol enables spatial redundancy to be introduced, which is an essential element for the reconstruction of a sky 



at a resolution greater than the focal plane resolution (without scanning) | Orieux et al. 2009 Champagnat et al. 
[2009] . 

2.1.2. Optics 

The Herschel Telescope is a classical Cassegrain instrument with a 3.5 m diameter primary mirror and a 308 mm 
diameter secondary mirror. The SPIRE photometer has three channels for a single field of view. The light is split 
by a combination of dichroics and flat folding mirrors. The spectral channels are defined by a sequence of metal 
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mesh filters and the reflection/transmission edges of the dichroics. They are centred at approximately 250, 350 and 
500 /im (noted as PSW, PMW and PLW respectively). We take the overall transmission curves of the wavelength 
filter /ife(A), for fc = 1, 2, 3, as given by the SPIRE Observers' Manual (no analytical form is available). 

The three detector arrays contain 139 (250 /im), 88 (350 /im) and 43 (500 /im) bolometers, each coupled to 
the telescope beam with hexagonally close-packed circular feedhoms. The beam solid angle is apodized by a 
bell-shaped weight whose width increases with A. Efforts have been made to correctly integrate the feedhoms in 
the instrument model but the detailed coupling of feedhoms on incoming radiation is, to the best of our knowl- 
edge I Griffin et al. 2002| , not fully understood at present. 



Our final choice as an effective PSF for the telescope coupled with feedhorns is a Gaussian shape ho{a, f3, A). 
This choice has two advantages: (i) it allows a closed equation for the instmment model (see Sec. |2.2| i, and (ii) it 
is in adequation with the global response measured from observations of Neptune (Griffin et al. 2010). As a first 
approach, we take isotropic Gaussians with standard deviations (To(A) = cA proportional to the wavelength since 
the width of the beam varies almost linearly with the wavelength. The coefficient of proportionality c is obtained 



by a least square fit of the Gaussian FWHM with the measured FWHM in [Ferlet 2007) . The widths obtained are 
close to the FWHM measured on the sky with 18.1", 25.2", and 36.9"at 250 /im, 350 /im and 500 /im, respectively 
(Griffin et al. 2010). The feedhorn diameter is 2FA, which introduces a focal plane sampling period of 2i^A (50" 
for the 350 nm array), or equivalently with sampling frequency fs ~ 0.02 arcsecond^^. 

The output after each feedhorn is then written as a 2D convolution of the input X{a, f3, A, t) and the effective 
PSF ho in addition to the wavelength filter 

A'<!"(A,t) = /ifc(A) [[ X{a,P,\t)K{a-ai^,P-Pi^,\)da<ip (3) 



where {aim, Pim) is the direction pointed by the feedhorn {l,m), for I = 1, . . . L and m = 1, . . . M. Finally, the 
optics is modelled as a linear invariant system w.rt. continuous variable. 

2.1.3. Bolometers 



To set up the bolometer model, we took the thermal model of | Sudiwala et al. 2002), which was also used in the 



simulator developed by | Sibthorpe et al. 2009) . Bolometers absorb all the received radiation 



P''"(i) = / X^"'{\,t)dX (4) 



A 



and this power provides the system excitation. The temperature T''™'{t) constitutes the system output to be deter- 
mined. The link between the input P{t) and the response T{t) is described by the differential equation deduced 
from a thermal balance. 



dT RiT)Vp Go . +1 _ 

where C is the heat capacity of the bolometer, R{T) is its resistivity, Tq is the temperature of the thermal bath, ly 
is a physical parameter that depends on the bolometer, Gq is the thermal conductance (at temperature Tq) and Vp 
and i?c are the polarization voltage and charge. No explicit solution of this equation is available in the literature. 
Sudiwala's approach | Sudiwala et al. 2002) , adopted here, is to linearize this equation around an operating point 



{T,P). In the following, we consider only the variable part of the flux and the constant part that defines the 
operating point is not included in the expressions. All the constants are defined with respect to the operating point. 
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For SPIRE, most of the observations should be carried out in the linear regime | Griffin 2006 2007 1. We thus 
consider that a first-order development is sufficient to model the bolometer behaviour correctly. Then, knowing 
how the resistivity R{T) varies with temperature, it is possible to determine the tension at the terminals. This 
first-order development models the bolometer as a first-order, low-pass linear invariant system having an impulse 
response 

= S'exp[-t/r] (5) 



where the gain S and the time constant r depend on the physical parameters in the differential equation | Sudiwala 



et al. 2002) . The values of these parameters are defined with respect to the operating point and correspond to the 



official SPIRE characteristics | |Griffin| [2006| |2007) . The output voltage around the operating point can then be 
written as a function of the incident flux: 

yim{t) ^ jJ^X'j^MKit' -t)dt' A\. (6) 

Finally, downstream, we have the read-out electronics, composed of several stages (modulation, amplification, 
low-pass filter, demodulation, quantification). However, except for the low-pass filters, they seem to have negli- 
gible effects relative to the other elements and are not taken into consideration in our model. The equations are 
nevertheless available | |Griffin| [2007) and it is possible to integrate them into the model. 

About the low-pass filters, they introduce a delay on the data with respect to the telescope position along the 
scan. As a trade-off between model accuracy and computation time, we have chosen to model the combination of 
the low-pass filter and the bolometer as a global first-order filter. The time constanj^ value (0.2 s) is taken to be 
representative of the combination. 

Finally we take into account regular time sampling that takes the values at times t — riTs (with a sampling 
frequency = l/T^ « 30 Hz) and then yj"^^ = yir,i{nT^), for n = 1, . . . TV. 

2.1.4. Complete equation of model 

Putting the above elements end to end gives the equation of the acquisition chain. For a spectral channel fc, the 
time signal at the bolometer (/, m) at time n is written 

^fc(A) / / X{a- pa{t).,l3 - Pi3{t),\)K{a - aim,P - A™, A) da d/3 h-b{nTs - i)dAdt . (7) 



This equation brings in four integrals: two coming from the optics (spatial convolution), one from the spectral 
integration and one from the time convolution. This is the fundamental equation of the instrument model since it 
describes the data yf^ bolometer by bolometer and at every instant as a function of the sky A" (a, /3, A). It should 
be noted that this model includes the discretization process in the sense that the data yf^ are functions of discrete 
variables (Z, to, n) and the sky A" is a function of continuous variables (a, /3, A). 



' Eventually, as for the illustration on real data Section 4.2 the correction of the low-pass filter can be performed using the 
Herschel Interactive Processing Environment (HIPE), and the time constant of the first-order low-pass filter is set to the time 
constant for the bolometer alone (5.7 ms). 
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2.1 .5. Model of sky for over-resolution 

The model of the sky is an important element for the reconstruction. As stated in the introduction and presented 



in Section 2.1.1 the spatial redundancy of the data should allow for partial reverse aliasing and enable a over- 



resolved sky to be estimated | Orieux et al. 2009 1 . The description of X must therefore be suitable for over- 
resolved reconstruction and, in particular, allow a fine description of the physical reality in connection with the 
above instrument model. 



Firstly, unlike conventional models | Sibthorpe et al. 2009 Cantalupo et al. 2010 1, we take into account the 
spectral variations of the sky within each channel. The emission observed by SPIRE is mainly due to dust particles 
in thermal equilibrium (between emission and absorption of UV and visible photons from incident radiation), and 
the intensities can be written 



^0 X 



(8) 



where T\g is the optical depth at wavelength Aq, (3 is the spectral index, B\ is the Planck function, and T the dust 
temperature. The SPIRE data alone do not allow the proper measurement of the dust temperature (the combination 



of SPIRE and PACS is mandatory, e.g., | Abergel et al. 2010| ), so we decide not to include the dust temperature in 
our sky model and work in the Rayleigh-Jeans approximation, so that B\{T) cx A^*. Moreover, we take (3 = 2, 
which is the "standard" value of the diffuse ISM (e.g., j jBoulanger et al.|ri996) ). Finally, we have 

X{a,f3,X) = X-'^Xia,P) . (9) 

with g = 6. However, as we will see in Section |2.2| the wavelength integration of the acquisition model will be 
done numerically. In other word, the spectrum profile is a degree of freedom that must be set in adequacy with the 
available knowledge about the observed sky. 

Secondly, X{a,P) is decomposed into a family of functions regularly shifted in space: ipijia^p) = tp{a — 
i6a, (3 — jSjs) where ip is an elementary function and {da, Sp) are the shifts between the in (a, (3). We then 
have: 

Xia,p)=Y,^vHa-^Sa,/3-jSp) (10) 

ij 

where is the decomposition function and x are the coefficients. 

One of the purposes of this decomposition is to describe the sky at a resolution possibly greater than the focal 
plane resolution. We therefore take an arbitrary nominal resolution corresponding to a target frequency width, 
greater than or equal to that of the data. To describe these maps, a natural approach is to choose tp as the cardinal 
sine for which the frequency width is the target width (and a step 6 equal to the inverse of this width) as described 



in Shannon's theory [Shannon & Weaver 1948 1. 

However, using cardinal sines requires analytical calculations that cannot be made explicit. To lighten the 
computational burden we will choose Gaussian ip functions from now on. These functions are parametrized by 
their spatial shifts {Sa, Sp) and their standard deviations (ctq., CT/j). The parameters (Jq, (5^3) are chosen to be equal 
to the inverse of the target width as for the cardinal sines. In the numerical treatments of Section]?] the values of 
{Sa, 6p) are equal to the pointing increment of the Large map protocol, i.e. 2" | Orieux 2009| . For the Gaussian 
function width parameters (ctq,, ct^), we determined the value that minimizes the difference between the width at 
half-height of the cardinal sine and the Gaussian: a^/p ~ 0.6 S^/p in a similar manner in a and /3. 
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2.2. Explicit calculation of acquisition model 

Given the linearity of the instrument model (|7]i and the decomposition ([9|)-(fT0l), the calculation of the model output 
for a fixed sky takes the form: 

Vhn^Yl^'i J ^^^f^kW jjj i}{a-i5a-Pa{t),li-j6fs-pi3{t)) 

ho{a - aim,P - Pim, A) dad/? K{nTs - t) dtdX . (11) 

Thus, to obtain the contribution of a sky coefficient Xij to a data item j/JJ^^, it is necessary to calculate four integrals, 
the discretization of which by brute force would result in very heavy numerical calculations. 

Concerning the optics, the convolution of the function ip with the optical response ho appears in Eq. ( [TT| i and, 
as these are Gaussians, the convolution can be written 



i>{a- iSa -pa[t),(i - jS/s - pi3{t))ho{a - aim, /? - Am, A) dad/3 



oc exp 



2S2 



(12) 



(Tq. This provides an exphcit resolution of the spatial 



with, in a similar manner in a and /3: S^y^ = ^a/p 
convolution. 

For the integral over time, only the constant velocity phases can be explicitly described for the Large map 
protocol. In order to integrate over time in ( [TT| ), we use the expressions of (|2]i for pa [t) and pp (t), which gives 



\-^hki\) I exp 



2S2 



exp 



{vpt + Cff +j6p - Pinif 

2E2 



K{nTs-t)<:\tA\. (13) 



It can be shown that explicit integration can be performed by including the Gaussians and the bolometer response 



(see details of the calculations in appendix [B] or in the supplements | Orieux 2010 1), and the model becomes 



S 



Vim — 



2V27rS 



-^^^ij / A '^ft.fc(A)exp 



(oq + nT^Va) {op + nT^vp) 



2S2 



2S2 



erfcx 



^\/2rE„ 

In this equation, the angles Oq, and op are defined by; Oq 



dA. (14) 



aim and op = cp + jSp - jSim- Moreover, 



y2 _ y2 2 
— ^p^a 



^a^p- 



The data point yf^ does not depend directly on the "scanning time" t as it is integrated and depends on time 
only through the sampling instants. These instants occur only in the form Oq + nTsVa and op + nTsVp, and the 
equivalent spatial samphng is thus finer than the equivalent angular separation of the horns. The introduction of 
this spatial over-samphng produces spatial redundancy, which, properly modelled and taken into account, is a key 
element for reconstruction that is over-resolved relative to the focal plane resolution. 

In addition, the time constant of the bolometer and the electronics t appears only in the term I]qS]^/(-\/2tE^), 
as an argument of the function erfcx. It is thus through this function erfcx that the influence of the bolometer and 
the electronics on the global spatial response of the system is taken into account. 
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The dependence on the wavelength through do (A) makes explicit integration with respect to A impossible. 
However, the integral depends neither on the data nor on the unknown object but only on the values defining the 
protocol. So, for a given protocol, these integrals can be calculated once and for all. Finally, the initial model 
brought in four integrations and the work described above makes three of them explicit. 

Eq. ( [I4] i models the acquisition of the data item y[Jjj at time n by bolometer {I, m) from the coefficients Xij 
that define the sky. These equations can be written 

y?r,i = X! '^Irnnii^tj) (15) 
ij 

where Ti. is calculated from Eq. ([14]). The model is linear and we can thus write 

y = Hx (16) 

where y and x are vectors of size LM N and IJ, and H is a matrix of LM N rows and IJ columns, each row of 
which can be deduced from ( [T4] i by varying I, m,n for fixed 



2.3. Invariant structure 

Initially, the physical model (|7]i is based on convolutive (so invariant) transforms w.rt. continuous variables. 
However, the discretization operation is irregular so the invariance property does not hold anymore. Nevertheless, 
the trace of this initial invariance continues to exist through the fact that if is a sum of terms at different spatial 



positions of the Gaussians (cf. Eq. ( 14 1). As the problem is now discretized, we seek to bring out an invariance by 
quantified shifts in 

+ riTsVa - aim 

for the a direction, and similarly for /?. 

Consequently, the bolometers are positioned on a hexagonal grid that can be decomposed into two rectangular 
sub-grids of step and Pg according to indices I and to. We thus have = IPa + toPq and the position of 
the bolometers during scanning becomes 

Oa + nTsVa ~ Ca+ iSa + nTsVa — IPa — mPp . 

If all the terms are multiples of a common factor Aq, the continuous shift is 

Oa + nTsVa — {no + ini + nn2 — In^ — mn4)Aa ■ 

This comes down to rounding the shifts to a multiple of A„. It can be said that the actual position of the bolometers 
is shifted or that the data are interpolated to the nearest neighbour The MADmap and SANEPIC methods use 
this idea but there is a notable difference: they perform the operation on a low-resolution grid, which results in 
limitation of the map resolution. In contrast, the developments proposed here exploit the idea of a high-resolution 
grid, enabling over-resolution or super-resolution reconstruction. By acting in the same way in the /3 direction, we 
have 

y?m — ^ij ^{i^o + + n-n2 — In^ — mn4)Aa, {n'^ + in\ + nn'2 — In^ — TOn4)A^ 
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and by computing the discrete convolution, we obtain 

y{t\f) = J2x,,H(^{t-t')A^,{j-f)A^') . (17) 

Thus y'll^ — y{i',j') if and only if 

i — i' = irii + In^ + m?T.4 — nn2 — Uq (18) 
j — j' = jn'i + In^ + mnn — nn'2 — Hq . (19) 

In these conditions, the data y, for a given scanning direction, are computed by discrete convolution ( [T7] l followed 
by (irregular) down-sampling defined by ([T8])-([T9]). 

This structure allows the model output to be computed much more efficiently. First of all, the decomposition 
by convolution then decimation is faster than the direct calculation and, what is more, the convolution can be 
calculated by FFT. Finally, given that only the impulse response is necessary, there is no need to calculate and store 
all the elements of the matrix. 

In this form, some calculations may be made even though they are useless, as the convolution is performed for 
all the indices whereas only some of them are used. The ratio depends on how fine the grid for the sky is relative 
to the size of the bolometer grid and the spatial shift. In practice, the excess calculation is reduced as we choose 
shifts {5a, Sfs) close to the delay between two sampling instants. Almost all the convolution results are observed. 

There is, however, the disadvantage that the bolometer positions are approximated. Yet these positions are 
important because they allow to get the best out of the data and to properly manage the information needed to 
estimate high frequencies. We choose a step A that is close to the minimum shift of the bolometers, Pq/15 at most 
for the smallest detector sampling step, i.e. about A w 2". The error introduced is thus small. This can be seen to 
be all the more valid when we consider the expected level of noise and telescope pointing errors, which are of the 
same order of magnitude, 2". 

Finally, the initial model ( [T6| l is decomposed in the discrete convolution defined by ( [T7| ) followed the (irregular) 
down-sampling defined by ([T8|-([T9|, that is to say H is factorized and: 

y = Hx = PH^x (20) 

where He is a convolution matrix and P a binary matrix that takes the values observed after the convolution. This 
is the pointing matrix similar to the one of SANEPIC and MADmap methods (but at a different resolution level). 
It has one, and only one, "1" per row since each data item can only come from one position. Some columns may 
be entirely zero as certain coefficients may not be observed. Conversely, some columns may contain several "1" 
because certain coefficients may be observed several times. From this, we also deduce that the matrix P is a sparse 
matrix. 

To sum up, by using an approximation on the position of the bolometers, we have separated the model H, 
which has no particular structure, into two sub-models H = PHc where He is invariant and P contains the 
non-invariant structure, which amounts to an irregular down-sampling. This decomposition is broadly similar to 
the one generally found in over-resolution in the field of image processing (see references in the Introduction). 

Fig. [1] presents this decomposition for the PSW detector with a velocity of 30 "/s towards the right: spatial 
redundancy contained in P (the blacker the pixel, the more it has been observed) and spatial impulse response (the 
time response of the bolometer and the electronics is clearly visible as the spatial extent of the Gaussian lobe). 
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Figure 1. Factorized physical model (PSW detector, velocity of 30 "/s towards the left): map of spatial redundan- 
cies P (left) and spatial impulse response He (right). The spatial scale are not proportional for better visualisation 
of the impulse response. 

2.4. Conclusion 

In conclusion, we have constructed a linear instrument model from the physical description of the phenomena 
involved during acquisition: scanning, optics, filters, bolometers and electronics are taken into account, together 
with a description of the sky in continuous variables in the three dimensions. We next explicitly described certain 
calculations and approached the model in a factorized form to lighten the numerical computation burden. 

The model presented here differs from those currently used in SANEPIC [Patanchon et al. 2008 1 or 
MADmap JCantalupo et al.| |2010) in that it takes the physics of acquisition into consideration. Moreover, un- 



like monochromatic models |Sibthorpe et al. 2009 1, the sky model extends spectrally over the whole channel. 



Again unlike jSibthorpe et al. 2009 1 , our bolometer model is linearized, which simplifies the developments and 
allows the bolometer time response to be made clearly explicit. 

Finally, the clear, consistent, global definition of the acquisition allows the spatial redundancy to be exploited 
directly and a processing method to be designed that uses these properties to estimate the sky at higher resolution 
than the focal plane resolution. 



3. Data inversion for liigli resoiution maps 

The previous Section was dedicated to the instrument model describing the relationship between the measured data 
y and the unknown sky X or its coefficients x through the equation 

y — T-LX + n = Hx + n . 

The matrix H is relatively complex and high-dimensional but the model ( fTSj l remains linear. The term n accounts 
for measurement and modelling errors and quantifies the data uncertainties. In this work, we consider that each 
bolometer denoted h is affected by an unknown offset Ofc. Eq. ([16]) can be rewritten for the bolometer h 

yt = HbX + Ob + nb. (21) 

where t/f, are data of bolometer b, Hi, is the corresponding part of the instrument model and Tif, accounts for eiTors 
at bolometer b. This Section presents the approach and methodology used to estimate the unknown x and the 
offsets o from the data y. 
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The question of sky reconstruction is a typical inverse problem and abundant literature is available on the 
subject fldier 2008 Demoment 1989 Tikhonov & Arsenin 1977 Twomey 1962| . As presented in the previous 
Section, the instrument model embeds convolutions and low-pass systems. Consequently, data are poor for reliable 
image reconstruction. The inverse problem is ill-posed and this is particularly true when over-resolution is intended. 
In this context, a naive inversion, such as a least squares solution, would lead to an unacceptably noisy and unstable 
solution. 



A usual class of solutions relies on regularization, i.e. the introduction of information | Idier 2008 Demoment 
1989 Tikhonov & Arsenin 1977 Twomey T962] to compensate for the lack of information in the data. A conse- 
quence of regularization is that reconstruction methods are specific to a class of map, according to the introduced 
information. From this standpoint, the present paper considers extended sources and relatively smooth maps. A sec- 
ondary consequence of ill-posedness and regularization is the need to balance the compromise between different 
sources of information. 



A possible solution relies on the TSVD | Hansen et al. 2000j that truncates the low singular values of H 
that are practically responsible for noise amplification. Another approach | |Richardson| |1972| stops the iterative 
algorithm before the solution is reached in order to prevent the occurrence of artefacts. The main drawbacks of 
these approaches are the unknown nature of the introduced information and the poor control (through a threshold 
or stopping criterion) of the trade-off with the data information. 

As already pointed out, a relatively spatially regular sky is expected. Since it is defined as a function of contin- 
uous variables, the regularity can be measured by the energy of derivatives of the function X. In the case of first 
derivatives in both directions, it can be shown (see appendix |A]) that 





2 


dX{a,l3) 


da 


+ 


dp 



(22) 



where D — + Dp \s obtained from the autocorrelation of the basis ij] and is similar to a discrete gradient 
operator. This relation illustrates the equivalence between the measure on the continuous function X and the 
measure on coefficient x, thanks to the use of Gaussian decomposition]^ 

With the regularity measure ( |22] i and the white Gaussian hypothesis for n, the solution is defined as the mini- 
mizer of the regularized least square criterion 



X — argmin Jx{X , o) 



where 



Jx{x,o) - \\v-nx~o\\ 



dX{a,/3) 


2 


dX{a,l3) 


1 


da 


+ 


d(3 





(23) 



(24) 



The parameter /i tunes the trade-off between adequation to the data and smoothness of the map. Since the set of 
ipij forms the basis of the function space, X can be expressed with a unique set x 



^ As an alternative, a non-quadratic norm of the derivative, e.g. convex penalties, could also be used. The interest of this is 
less penalization of high gradients in the map. Unfortunately, an explicit measure on coefficients, as in the quadratic norm, is 
not explicit. 
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The criterion Eq. ( |24] i is defined on X but the same solution is obtained with an equivalent criterion on the coeffi- 
cient 

x,o = arg min {x^ o) 

where 

Jx{x,o) ^\\y - Hx - oW"^ + ^lx^Dx (25) 

and the property 

J^{x,d) = Jx{X,d). (26) 



Remark 1. - A Bayesian interpretation of criterion ( p5| l is the Gaussian posterior law 

p{x,o\y) (xp{y\x,o)p{x)p{o) 
oc exp 



.^\\y-Hx-or-^x'Dx 



with Gaussian iid likelihood, Gaussian correlated prior and flat prior law on E, for o. Consequently, the minimum 
of the criterion is the maximum of this law 

x,o = arg maxp(a;, o\y) 

x,o 

where /i — a'^/a'^. An advantage of the Bayesian interpretation is the ability to derive an uncertainty around the 
maximum through the variance (see Sec. |4| of the posterior law. Another important advantage of the Bayesian 
interpretation deals with the estimation of p and of instrument parameters and we have a second paper devoted to 
this point | Orieux et al. 2010) . 



The proposed algorithm for the computation of x and o is an alternate minimization algorithm: after an initial- 
ization, the following two steps are repeated 

1 . Find X for fixed o 

x'' = argmin \\y - Hx - o'^jp + fix^Dx (27) 

X 

2. Find o for fixed x 

o'=+i = arg min \\y - Hx'' - o| | ^ (28) 

o 

until a criterion is met. For fixed x, the solution is straightforward and Ob is the empirical mean of the residual 
j/b — HjjX for each bolometer separately. For fixed o, the solution Eq. ( p7] i is unique (because the criterion is 
convex) and explicit 

X = {WH + ^iD)~'^ H\y- o). (29) 

The estimator is linear w.rt. data y. Unfortunately, since H is not circulant, x cannot be computed with a "brute 
force" algorithm: the practical inversion of the Hessian matrix H^H + fxD is impossible (the size of this matrix 
is the square of the number of coefficients x). The proposed solution relies on an iterative conjugated gradient 
descent algorithm [Nocedal & Wright] |2000[|Shewchukl|1994[ . Because of the convexity and differentiability of 
the criterion, the algorithm is guaranteed to converge to the solution. 
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The conjugated gradient is a first order algorithm that operates the gradient gk at current point Xk 

Qk = 2H\y - o) - 2{H'H + ^lD)xu (30) 



to compute a new point. As illustrated by Eq. ( |30| l, the most expensive part is the computation of the product 
between the matrix H*H and the current point X}.. 

Remark 2. - As described in appendix |Aj the regularization part Dxk is computed by FFT. Moreover, as de- 
scribed in appendix |c] the model-related part H^Hxk can also be efficiently computed based on FFT, decimation 
and zero-padding. 

4. Experimental results 

This part illustrates the improvement that our approach can bring, first using simulated data and then with actual 
data transmitted by Herschel. 

4. 1. Simulated data 

4.1.1. Experimental protocol 



We chose three 20' x 20' maps used by the SPIRE consortium to assess reconstruction methods | Clements et al. 



2006): (i) a map of galactic cirrus (Fig.|4]i complying with the a priori regularity model, (ii) a map of galactic cirrus 



superimposed on point sources (Fig.|7]i and (iii) a galaxy map (Fig. [8]). 

The study concerns the PMW channel and the Large Map protocol with three scans in each dkection and a 
velocity of 30"/s. The data were generated by the instrument model ( fTSj l, considering for this simulation part that 
the bolometers are not affected by any offset. Moreover, the sky spectrum profile is a degree of freedom of our 
acquisition model since the wavelength integration is numerical (Eq. (|T4)l). The profile have to be set in adequacy 



with available knowledge about the sky, see Section 2.1.5 We assume for the simulations and the inversions that 



the sky spectrum is flat (g ~ 0). The noise is zero-mean white Gaussian noise and we consider three levels 
characterized by their standard deviation cr„ ("standard noise" hereafter), 10cr„ ("high noise") and 0.1 cr„ ("low 
noise"). The standard deviation is the same for aU the bolometers and, unless stated otherwise, all the data sets 
were generated with the same noise realization. 

The proposed reconstruction for the 20' x 20' maps are performed using 6a = = 2", i.e. maps of 600 x 600 
coefficients. Our results are compared with the map obtained by coaddition, with 6" as pixel size. 

In accordance with what was presented in Section[3] the map is reconstructed as the minimizer of criterion ( |23] l- 
(|25| and the minimization is performed by a conjugate gradient algorithm with optimal step size. Fig.|2]illustrates 
the behaviour of the algorithm for cirrus in the "standard noise" case. As expected, the value of the criterion 
decreases at each iteration and stabilizes. A few tens of iterations appear to be sufficient to reach the minimum in 
the example presented and in the other tested examples. 

In the simulated cases, as the original map (the "sky truth") is known, we can make a quantitative assessment 
of the reconstruction through an error measure defined by: 

e = El4-^l/ EI4I (31) 
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Figure 2. Value of the criterion during iterations of the minimization algorithm (conjugate gradient with optimal 
step). Case of cutus with "standard noise". 



where a;* is the true map and x is the reconstructed map. These measures enable the eiTors produced by the 
different methods to be compared quantitatively. 




io'° 10" 



Figure 3. Reconstruction eiTor e vs regularization parameter /i. Case of ciiTus with "standard noise". 



In the case of the proposed reconstruction, the estimate depends on the regularization parameter: x — x{ij.). 
Fig. |3]illustrates this dependence through the error e for cirrus in the "standard noise" case. A non-zero optimum 
value /iopt appears (here ^ 10^^) for which e is a minimum, thus confirming the interest of the regularization. 
A value lower than 10^^ leads to an under-regularized map and a value greater than 10^"^ to an over-regularized 
one. In the following, it is, of course, the optimal value that is used to reconstruct the maps. Also, as far as the 
map sensitivity to the regularization parameter is concerned, it appears empirically that /i needs to vary by a factor 
2 around fiopt for a visible modification of the map to be obtained. This result is confirmed in Fig. |3] where the 
minimum is not very marked relative to the horizontal scale. Fig. [3] also quantifies the improvement provided by 
the regularization: the errors for the non-regularized and optimal-regularized maps are 0.12 and 0.08 respectively, 
i.e. a gain of 33%. 



4.1 .2. Restoration of galactic cirrus 

Fig.[4]sums up the first results concerning the cirrus in the "standard noise" case. The proposed map is very close 
to the true one. In particular, our method restores details of small spatial scales (with spectral extension from low 



to high frequency) that are invisible on the coaddition but present on the true map (see also the profiles in Figs. 4(d) 



and 4(e) 1. In particular, the fluctuations around pixels 250 and 350 are well restored. In addition, our method also 
correctly restores large-scale structures, corresponding to low-frequencies down to the null frequency (mean level 
of the map). We conclude that our method properly estimate the correct photometry. 
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Figure 4. Comparison of results. Fig. 4(a) shows the true map. Fig. 4(b) presents the proposed map and Fig. 4(c) 



the coaddition. A horizontal profile is shown in Fig. 4(d) and Fig. 4(e) gives a zoom. 



Remark 3. - Moreover, the reconstruction method is linear with respect to the data (see Section|2]), which means 
that the use of arbitrary units is valid. 




10"^ 10"' 10"^ 10"' 10"^ 10"' 
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(a) Standard noise (b) High noise (c) Low noise 

Figure 5. Circular means of power spectra for the three levels of noise (standard deviations: cr„, 10 (Jn and 0. 1 (t„). 



To quantify the gain in correctly restored bandwidth, we look at the power spectra of the maps (Fig. [5]l for the 



true sky, the sky convolved with the PSF, the coaddition and the proposed sky. As mentionned in Section |2.1.2| 
the sampling frequency of the focal plane is /g ss 0.02 arcsecond^^. Consequently, by the Shannon theorem, the 
acquired data during one integration cannot correctly represent frequencies above fs/2 « 0.01 arcsecond^^. We 



have also seen in Section 2.1.2 that the FWHM of the PSF is 25.2"at 350 /im, i.e. a cutoff frequency of the optical 
transfer function of « 0.04arcsecond^^. The attenuation effect of the convolution by the PSF on the true map 
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(or of the multiplication of the Fourier transform by the transfer function) is visible on the power spectra of the 
convolved and coaddition maps, for all frequencies above w 0.008 arcsecond^^ (Fig.jsj). 

Regarding our approach, the power spectra of the proposed map perfectly follows the power spectra of the 
true map, from the null frequency up to a limit frequency that depends on the noise level. In the "standard noise" 



case. Fig. 5(a) this limit is 0.03 arcsecond"^ that is to say three times the limit frequency of each integration 
(/s/2 « 0.01 arcsecond^^). It illustrates that our method also takes full advantage of the high frequency temporal 
sampling. In any case and compared to the coaddition, we have multiplied by a factor w 4 the spectral bandwidth 
(starting from the null frequency) where frequencies attenuated by the optical transfer function are accurately 
corrected. The high-resolution reconstruction potentiality of the proposed method is clearly demonstrated in this 
experiment. 
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Figure 6. Uncertainty provided by the a posteriori standard deviation a. Fig. 6(a) shows the map of the standard 



deviation for each pixel and Fig. |6(b)| gives a profile. Fig. |6(c)| shows a profile of the true map as a solid line and 



the two dashed line give a ±3(t interval around the estimated map. Fig. 6(d) shows the PSD of the true map as a 
soUd red line and the two dashed hne give a ±a interval around the estimated PSD map in the "standard noise" 
case. 



Our method also gives a value for the uncertainty on each pixel through the Bayesian interpretation presented 
in Remark [T| This uncertainty is provided by the standard deviation a of the a posteriori law and a map of these 



standard deviations is shown in Fig. 6(a) The uncertainty is naturally smaller in the centre of the map and grows 
as we move away from the centre because the data contain less information. Fig. |6(b)| shows a profile of the 



uncertainty and Fig. 6(c) gives a profile of the true map and a ±3ct interval around the estimated map. We thus 
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show that the true map is indeed inside the interval around the estimated one. The Fig. 6(d) gives the true PSD and 
a ±1(7 interval the estimated PSD. Again, up to the 0.03 frequency limit, the true PSD is inside the interval. 

The possibilities of restoring frequencies obviously depend on the noise levels. Fig. |5] compares the spectra 
obtained with the three noise levels, the parameter /i being chosen to be optimal each time from the point of view 
of the error Eq. ([31). When the noise level is lower, it is possible to restore slightly higher frequencies: up to 
0.03 arcsecond^^ for "low noise", as against 0.025 arcsecond^^ for "standard noise". Conversely, in the case of 



"high noise", our method no longer restores the frequencies attenuated by the optical transfer function Fig. 5(b) 
The deconvolution effect is reduced and the essential effect is one of denoising. Nevertheless, the proposed method 
gives better (or equivalent) results than coaddition in all cases. 

4.1 .3. Other types of sky 




Our method is based on spatial regularity information. To assess how robust it is, we tested it with two other 
types of sky in which the spatial regularity is less pronounced: galactic cirrus superimposed on point sources, and 
a galaxy image. The results obtained are presented in Figs. [7] and [8] 



For the first case, the characteristics of the reconstructions can be seen on the maps by comparing Fig. 7(b) and 



Fig. 7(c) the coaddition map is smoother than the one we propose. In addition, some of the point sources of the true 
map are visible on the proposed map but not on the coaddition one. If we observe the profiles presented in Fig. |7(d)| 
we note that the point source near pixel 350 is underestimated but markedly less so by the proposed method than 



F. Orieux et al.: Over-resolution: instrument model and regularized inversion. 



19 



,10"- 



,10" 



(a) True map 
-xlOl' 



(b) Proposed 













— True 












— Proposed 
— Co-add 
















) 





























,10"- 



I 




(c) Coaddition 



-200 200 400 600 
(d) x-th profile 




TOO 150 200 250 300 350 
(e) x-th profile (zoom) 



I 



Figure 8. Restoration of galaxy. 



by coaddition. Rebounds also appear around the point sources, a feature characteristic of linear deconvolution 
(resulting from a quadratic criterion) in the presence of point sources. 

The galaxy does not contain point sources but has spatial structures that are more complex than the galactic 
cirrus. These structures are considerably better restored by our method than by coaddition and the comparison 
between Fig. |8(b)| and Fig. |8(c)| is instructive from this point of view. In particular, the double structures in the 



arms of the true galaxy (Fig. 8(a) i are correctly separated by the proposed method (Fig. 8(b) i but not by coaddition 



(Fig. 8(c) I. Moreover, while the profile of Fig. 8(d) shows real restoration of the fluctuations between pixels 200 
and 400 by both the proposed method and coaddition, the proposed restoration is clearly superior to coaddition. 

In conclusion, the proposed method is relatively flexible and shows a good restoration capacity for various 
types of map. In particular, it possesses a certain robustness with respect to an input sky presenting characteristics 
that are poorly taken into account by the a priori model based on regularity information. It provides a sky that is 
closer to the real one than that obtained by coaddition, even in the least favourable cases. 



4.2. Processing real data 

We conducted the first tests with real data of the reflection nebula NGC 7023 and of the Polaris flare (which is a high 
Galactic latitude cirrus cloud) performed during the science demonstration phase of Herschel and already presented 
in I Abergel et al. 2010 1 and |Miville-Deschenes et al. 2010| , respectively. In order to run our algorithm, we took 
the Level- 1 files processed using HIPE. The true sky is not known, of course, so the value of the regularization 
parameter was fixed for each of the spectral channel by a trade-off between the gain in spectral resolution and the 
amplification of the 1// noise. 
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Figure 9. Central part (23'x23') of NGC7023 in the three channels. Top panels: coadded maps; bottom panels: 
proposed map. 
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Figure 10. Brightness profiles along the three sections shown in the top left panel of Fig. |9j Each panel shows 
the profiles along one section within the PSW, PMW and PLW channel, offset for clarity from bottom to top 
respectively. Left panel: horizontal profile; central and light panels: vertical profiles. Coefficient are 2"spaced. 
Black: coadded maps, blue: proposed maps. 

Figs.[9]to[T2]illustrate the results for NGC 7023 and the polaris flare. The gain in spatial resolution is spectacular 
in the three channels. It is interesting to note that the map of NGC 7023 obtained by our method in the longest 
wavelength channel (500 /im, the PLW channel) shows spatial structures that are not visible in the coaddition but 
are real since they are visible at shorter wavelengths, as illustrated for instance on the last panel of Fig. [TO] The 
same panel also shows that negative rebounds appear on the sharpest side of the brightest filament of NGC 7023. 
This filament is the naiTowest structure of the map and its width is comparable to the width of a point source. 
Comparable rebounds were also seen on our simulations with point sources (Fig. [7]). The Polaris flare does not 
contain comparable bright and narrow filament, so the proposed map does not present this kind of artifact. 
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Figure 11. 50'x50'square in the Polaris flare (the total field observed during the science demonstration phase of 
Herschel is presented in | |Miville-Deschenes et al.|[2010[ ). Top panels and from left to right: coadded maps in the 
PLW, PMW and PSW channels, respectively; bottom panels: proposed results in the three channels. 



Fig. 12 shows that the power spectrum of the proposed map of the Polaris flare in the PSW channel follows 
the expected power law typical of the infrared emission of high Galactic cirrus P{k) a k'' with 7 = —2.7 (e.g., 
I Miville-Deschenes et al. 2010| ) on a frequency range from 10^"^ arcsecond^^ to 3 x 10^^ arcsecond^^. As for 
simulated data. Section [4T| the attenuation by the Optical Transfer Function (OTF) is accurately corrected up to 
the frequency where the noise is dominant. Thanks to this correction, the contrast of the small scale structures in 



the maps is enhanced (Figs. 1 1 and 12 1, since the energy of each structure is put in a smaller number of pixels than 
for the co-added case. 



5. Conclusion 

We have proposed a new generic method of over-resolved (with respect to the detector resolution) reconstruction 
of images for scanning instruments and have applied it to the SPIRE instrument of the Herschel observatory. 

The first key element is an instrument model that precisely describes the physical processes involved in the 
acquisition. To explain the data in a reliable way, our model combines the descriptions of three elements: (i) the 
sky as a function of continuous variables in the three dimensions (two spatial and one spectral), (ii) the optics 
and the rest of the instrumentation (bolometer, electronics,. . . ) and (iii) the scanning strategy. We thus arrive at 
a linear model in integral form (Eq. (|7]i). We then write this as a matrix expression (Eq. ( [T6| )) by making certain 
calculations exphcit. Next, by coming close to the pointed positions (on a high-resolution grid), we decompose it 
into a convolution followed by (irregular) down-sampling (Eq. (|20)i). This model provides the most faithful link 
between the data, the sky actually observed, and the instrument effects. 
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Figure 12. Circular means of the power spectrum of the 50' x 50' square shown on Fig {Tl]and in the PSW channel. 
The red line shows the power law P{k) oc fc^ adjusted on the data in a frequency range from 10^'^ arcsecond"^ 



to 3 X 10 ^ arcsccond ^, with 7 = —2.7. At smaller frequencies, | Miville-Deschenes et al. 2010 1 have shown that 



the SPIRE spectra are attenuated compared to IRAS, which is likely due to the correction of 1// noise attributed 
to thermal drifts in the preprocessing of the data. The pink soUd Une shows the Optical Transfer Function (OTF). 



On the sole basis of this instrument model and the data, the inversion is an ill-posed problem. The lack of infor- 
mation brought by the data, considering the limitations of the instrument, leads to instability of the inversion, which 
is all the more noticeable when the target resolution is high. This difficulty is overcome by a standard regularization 
method that constitutes the second key element. The method relies on spatial regularity information introduced by 
quadratic penalization and on a quadratic data attachment term, the trade-off being governed by a regularization pa- 
rameter. Thus, the inversion is based on a relatively standard linear approach and its implementation uses standard 
numerical optimization tools (conjugate gradient with optimal step). 

The presented results for the SPIRE instrument give a spectacular illustration, on simulated and real data, of 
the potential of our method. Through the use of the accurate instrument model and a priori regularity information, 
we restore spatial frequencies over a bandwidth ^ 4 times that obtained with coaddition. In all channels, the 
attenuation by the optical transfer function is accurately corrected up to the frequency where the noise is dominant. 
The photometry is also properly restored. 

There are two main straightforward perspectives to our method: firstly, the instrument parameters must be 
known and, secondly, the regularization parameter must be arbitrarily chosen. Future work will propose a "non- 
supervised" and "myopic" inversion method for setting these parameters automatically. 

As other perspective, quadratic prior adequacy is known for possible excessive sharp edges penalization in the 
restored object. The use of convex L2 — Li penalization |Kiinsch 1994[ Charbonnier et al. 1997) can overcome 
this limitation. Secondly, estimation of parameters of correlation matrix (cutoff frequency, attenuation coefficients, 
spectral profile, . . . ) could be achieved for the correlation matrix of the object or the noise (typically for the 1 // 
noise). 
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Finally, a really valuable perspective could be to introduce the spectral dependence between the different chan- 
nels in the data inversion. The conjunction with a PACS direct model and the joint inversion of SPIRE and PACS 
data would greatly improve the map reconstruction. 
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Appendix A: Energy spectral density 

This appendix gives the details of the calculations concerning the regularity measure and its frequency interpreta- 
tion, which are used in Section [3] Based on the decomposition of Eq. ( fTO] ), the energy of the first derivative can be 
written 



dX 
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da 


J 





dad/3 



iji'j' 
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da 



A'j' 



d_ 
da 



) dad/3 



By noting the derivative i/j'^ = dijj/da, we bring out the autocorrelation = * of the first derivative of 
the decomposition function and we have 
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As there is a finite number of coefficients xij, the measure can be put in the form of a quadratic norm 

2 

— - fJC D {y 3C 



(A.l) 



dX 



da 

where the matrix Da is obtained from ^Pq. Considering the invariant structure of ( |A.l| i, the matrix Da has a 
Toeplitz structure. The calculation is performed by discrete convolution and can be done by FFT. 
By introducing the dimension /?: 

dX ^ dX ^ . . 
— + — =xDaX + xD,x 

The quadratic regularity measure on the function X with continuous variables is expressed through a quadratic 
regularity measure on the coefficients x. 

The autocorrelation Fourier transform (FT) is the energy spectral density, i.e. the squared modulus of the FT 



*a(/a,//3) 
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i;'a{a,P)e-'^^''^°'f-+'^M ^ad(3 
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where ^ is the FT of ijj. When the dimension j3 is introduced, the a priori energy spectral density for the sky 
naturally has circular symmetry 

"^{fcJp) = 4^' {fl + /|) |V'(/a,//3)f . (A.2) 

This calculation brings out the frequency structure introduced a priori for the sky according to the chosen function 
■i/j. This is a high-pass structure since the factor + tends to cancel ^ around zero, which is consistent with a 
regularity measure. 



Appendix B: Explicit calculation of the model 

In order to integrate over time in ( [T3| ), we use the expressions of (|2]i for (t) and (t), which give 
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With the bolometer response 
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with Oq, = Cq, + a'-' — a;m and 0/3 = cp+ (3^^ — film ■ This is the integration of a truncated Gaussian as the argument 
of the exponential has a form that is quadratic in t. 



B. 1. Calculation of the argument 

Here, we express the quadratic form in question 

iTEjiva^t + Oo,)^ + T^l{vpt + opy ~2^l^t _ 1 n{t) 



2WIt 



Developing and factorizing the numerator n{t) gives it the form 

{t + af + h- 



n{t) = tY^I {vie + 2vc,0gt + ol) + tI]2 (vjf + 2vpOpt + o^) - 2Y.mt 



S2 



with the constants a = Y? {rY.'j^VaOa + TYlv^op - and h = rS^ (^S^^s ^ Y^oj^ and also 

T {ri^pvl + Yilvj^ . Putting this i-quadratic form into the integral we obtain 
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where the function erf is defined by: 
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erf (a;) = — / e"®' de* = -erf (-x) 
V^r Jo 

This expression can be simplified by using the function erfcx(a;) — exp(a;2)(l — erf (x)). 
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B.2. Argument of the exponential 
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To make the written form lighter, we can set S = E^S^S|. The factor intervening in the function exp is 
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So, by injecting this expression in (|B.2|i, the function erfcx appears 
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The values of 5, a and b can be replaced. First of all, the argument of the exponential is written 
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(B.3) 



and the terms nTs/r simplify. We then use the expression for 

„2j.2 _ n^T^r [^vl + J:lvfj _ „2y2„2 n^r^ 
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to bring out two perfect squares. Finally the argument of the exponential ( |B.3| l in ( |B.2| l is written 

(oq + nTsVaf {of} + nTsVpf 



2E2 



2S^ 



(B.4) 



which is exactly the argument of a bivariate Gaussian. We again find the same standard deviations and E^. 
However, the response of the optics, initially (oa, o^), is now shifted by {nTsVajnTsVp), i.e. the pointing difference 
between two successive time samples. 



B.3. Argument of the function erfcx and final expression 

Another term is needed if we are to know the global response. It comes from the function erfcx, which corresponds 
to the influence of the bolometer. The argument of the function erfcx is 

■^2'sn 

/3 



nTs +a TiTs + EM T^VaOa + rT^i^vpOp - E^E 
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(B.5) 



where E^ = E|u^ + E^i;^ and what is of interest here is that the same factors are found in the argument of the 
exponential. To know the global response, we need to bring everything together. By injecting the expressions of 
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the arguments ( |B.4| l and ( |B.5| l, we obtain 
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(B.6) 



with, similarly for a and /?: ,„ = '^Ii/b + '''o' which finishes the integration of ( [T3| ) over time. 



Appendix C: Direct modei computation aigoritiim 



This part gives some more details on the concrete calculation of a model output Hx of Section 2.3 First of all, 
there are four different impulse responses whatever the number of scans. For two scans in the same direction, the 
convolution is the same. Thus we can construct four different convolution matrices Hi for i = 1, 2, 3, 4 and apply 
four different discrete convolutions to the coefficients x. 

We can also deduce the structure of the transpose of the model ~ H^P*^. The matrix P* is a data sum- 
mation / zero fill matrix (addition of the data that possess the same pointing while setting the other coefficients to 
zero), and ff* corresponds to a convolution with the space reversal impulse responses. 

The product by P* is very similar to the construction of a naive map except that the data are added rather 
than averaged. Also, the operation is done by velocity and not globally. Finally, the products by He and iJ* are 
convolutions implemented by FFT 
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